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We develop the formalism for determining the quasinormal modes of general relativistic multi- 
fluid compact stars in such a way that the impact of superfluid gap data can be assessed. Our 
results represent the first attempt to study true multi-layer dynamics, an important step towards 
| considering realistic superfluid/superconducting compact stars. We combine a relativistic model 

. for entrainment with model equations of state that explicity incorporate the symmetry energy. Our 

analysis emphasises the many different parameters that are required for this kind of modelling, and 
the fact that standard tabulated equations of state are grossly incomplete in this respect. To make 
progress, future equations of state need to provide the energy density as a function of the various 
nucleon number densities, the temperature (i.e. entropy), and the entrainment among the various 
' components. 
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I. INTRODUCTION 



The fluid approximation [l[ is a necessity for modeling systems containing so many elementary constituents that 
(on a macroscopic scale) they form a continuum, and can redistribute energy and momentum among themselves. The 
model is based on the notion of fluid clement, which is small enough to be infinitesimal with respect to the system 
en masse, but large enough to contain, say, an Avogadro's number worth of particles. In a superfluid system one 
must consider different, dynamically decoupled, yet co-existing fluids. Each individual fluid has its own collection 
of fluid elements and each spacetime point in the system will have as many fluid element worldlines passing through 
it as there are independent fluids. The extent to which different fluids are coupled depends largely on the available 
^— «) | dissipation mechanisms, eg. friction due to interparticle scattering. In a multi-fluid system this is a complex problem 
0. Generally, the system dissipation depends on how energy and momentum are exchanged as the fluid elements 
expand and contract, how they slide across each other, how they rotate about each other, and how they flow through 
each other. 

Neutron stars are believed to be prime examples of general relativistic, multi-fluid objects. To date, on the order of 
a couple of thousand pulsars have been observed [|[. Yet, an extrapolation of the local population data for our galaxy 
suggests the existence of about 1.6 x 10 5 normal pulsars, around 4 x 10 4 millisecond pulsars, and about 2 x 10 8 neutron 
stars that are no longer active pulsars. Most of these should be extremely cold in the sense that their temperature is 
several orders of magnitude less than the Fermi temperatures of the independent (massive) particle species, i.e. 10 12 K. 
In fact, neutron stars cool to temperatures below 10 9 K relatively soon after birth. This is the expected transition 
temperature 0, [H, Q for neutrons and protons to become superfluid and superconducting, respectively. We therefore 
anticipate that most neutron stars in our galaxy have at least two superfluids/supcrconductors in their cores. In fact, 
neutron superfluidity is a key ingredient in most models of large pulsar glitches @, 3 , with (catastrophic) transfer of 
angular momentum via vortices from a superfluid component to the crust leading to the observed spin-up. 

Ever since neutron star superfluidity was first suggested [9|. we have seen a concerted effort aimed at developing 
our understanding of the various phases of matter [l(| El, E3 • This has led to a picture where a typical neutron 
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star has a number of distinct "layers". From the outer to the inner crust, protons are locked inside increasingly 
neutron-rich nuclei which are embedded in a degenerate normal fluid of electrons. At the base of the crust, the nuclei 
are believed to assume exotic "pasta" shapes [13|. Moreover, in the crust region the long-range, attractive component 
of the nuclear force should lead to Cooper pairing of neutrons in x So states. The crust nuclei will thus be embedded 
in and penetrated by superfluid neutrons. Short-range repulsion in the nuclear force and the spin-orbit interaction 
allow neutrons to pair as 3 P2 states in the more dense regions of the outer core [HI]. There are no nuclei in this 
region and the protons that remain are dilute enough to feel only the long-range attractive part of the nuclear force. 
Hence, they are expected to pair in states. The core neutrons and protons are embedded in a highly degenerate 
normal fluid of electrons. At high enough density it becomes energetically favorable for more massive particles to form 
(eg. muons in lieu of electrons). At the most extreme densities, quarks may become deconfmed, possibly opening up a 
number of channels for attractive interactions and hence many Cooper pairing possibilities (e.g. so-called CFL matter 
, Gil). Each individual Cooper pairing will lead to individual condensates, and potentially many inter-penetrating 
fluids. The actual number of dynamically distinct fluids depends very much on the details of dissipation @, [l3|- The 
scattering time-scale must typically be much greater than the characteristic dynamical time-scales in order for fluid 
components to decouple. 

A moderately realistic neutron star model must account for the presence of different regions in the star. Describing 
a neutron star as a multi-layer system, one should at the very least account for the presence of the (multi-fluid) 
superfluid/superconducting regions and the elastic crust component. Even though there has been progress in this area, 
we are still not at this level of sophistication. The current state-of-the-art is represented by (llj , where quasinormal 
modes for a core-envelope model are calculated, and a recent study of axial crust oscillations in the relativistic Cowling 
approximation [l9j . The aim of the present work is to improve on the first aspect of the modeling. We will include, 
for the first time, the physics of the superfluid phase transition in the construction of the spherically symmetric and 
static background model. The motivation for this is that the different regions of superfluidity may not overlap [|| , 
and the fluid dynamics will change depending on whether a given layer is normal or superfluid. This represents 
an important improvement on previous models, since it allows us to quantify how changes in the superfluid energy 
gap affect the quasinormal-mode spectrum. We leave inclusion of the crust elasticit y fo r future studies. A fully 
relativistic formalism for a crust penetrated by a superfluid has been developed, see eg. [20(, but it is yet to be applied 
to neutron star dynamics. 

In a mixture of the supcrfluids He 3 and He , it is known that a momentum induced in one of the constituents will 
cause some of the mass of the other to be carried along, or entrained pll . [22j ■ Another example is the entropy (which 
can be considered as a massless fluid). In fact, in superfluid He 4 the so-called "normal" fluid density is directly 
proportional to the entrainmcnt between the atoms and the entropy. In a neutron star, the strong interaction leads 
to entrainment between neutrons and protons. Entrainmcnt is a multi- fluid effect, that has no counterpart in a 
single- fluid system. When the neutrons start to flow they will, through entrainmcnt, induce a momentum in the 
protons and subsequently the electrons (23j . In principle, there could be entrainment between each and every fluid of 
a multi-fluid system. In this analysis we will only consider the entrainment between neutrons and protons, for which 
we use the fully relativistic mean field model developed by Comer and Joynt [24|. 

The outline of this paper is as follows: In Sec. [IT] we give details on the two-fluid formalism, and the set of equations 
used to model the quasinormal modes. Sec. IIIII discusses an expansion for the local matter content that is adapted to 
include entrainmcnt at the appropriate order. Sec. IIVI gives the specifics on the local matter content, i.e. the equation 
of state, the gap data, and the relativistic entrainment. The following Sec. IVl provides the results of our analysis of 
modes for two equations of state, with or without entrainment, and different "temperatures" of the star. Sec. IVII 
gives some concluding remarks and discusses to what extent the extant literature and computational infrastructure for 
compact object equations of state are adequate for supporting our kind of analysis. Finally, the Appendix describes 
the technique used to obtain the results in Sec. [Vj We use "MTW" conventions throughout. 

II. GENERAL RELATIVISTIC TWO-FLUID FORMALISM 
A. The full formalism 

We will use the formalism developed by Carter, Langlois, and their various collaborators [25|, (2(| [27], [H, HI, HO, 
HH, [H, [H, (see Andersson and Comer [l[ for a recent review) . The fundamental fluid variables consist of the 
conserved nucleon density four-currents, to be denoted nf^ where x = {n, p} is a so-called constituent index (which is 
not summed over when repeated). From the currents can be formed three scalars: = —g^n^n^, = —g^n^n^, 
and n^ p = — g^n^n^. Given a master function — A(n^, np, n^ p ) (the two-fluid analog of the equation of state), the 
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stress-energy tensor is 



where 



is the generalized pressure and 



* = A - - ng^S 



= g M „ (£ x < + A* y n v y ) 



(1) 



(2) 



(3) 



is the chemical potential covector. It is also the momentum canonically conjugate to the current n£. Formally, the 
A xy and B x coefficients are obtained from the master function via the partial derivatives 



A xy = A yx 



dA 

dnl y 



B x 



; dA 

dnl 



(4) 



The fact that the momentum fj, x is not simply proportional to the corresponding number density current n£ is a result 
of entrainment. It vanishes if the A xy coefficient is zero. 

Finally, the equations for the neutrons and protons each consist of a conservation equation 







and an Eulcr equation 



where the vorticity two-form is defined by 



(5) 



(6) 



(7) 



The square brackets indicate antisymmetrization of the enclosed indices. Comer |33| and Prix et al [35| discuss in some 
detail why the assumption of separate conservation laws for the two fluids should be reasonable for slow rotation and 
quasinormal-mode calculations (excluding the thin transition layers that separate single fluid and multi- fluid layers). 
Note that the above way of writing the Euler equations makes manifest its geo metric meaning as an intcgrability 
condition for the vorticity, a point that has been much emphasized by Carter [25l | (see also [H). 



B. Equilibrium models 

In order to determine the background fluid configuration we need to evaluate the associated metric. We take our 
equilibrium configurations to be static and spherically symmetric. The metric can thus be written in the Schwarzschild 
form 

ds 2 = -e v dt 2 + e x dr 2 + r 2 (d9 2 + sin 2 0# 2 ) . (8) 

The required metric coefficients are determined from two components of the Einstein equations, which can be written 

1 — e A 1 — e A 

A' = — - 87rre A A , v = + 87rre A * . (9) 

r r 

A prime represents a radial derivative, and it is to be understood that A = A(n 2 ,n 2 ) and \& = ^(n 2 ,?! 2 ) in the 
two- fluid layers and A = A(n 2 ) and \& = ^(n 2 ), where n = n n + n p is the total baryon number density, in the 
single-fluid layers. 

The equation that determines the radial profile of n x (r) in the supcrfluid layer is (30j 



where 



B*°<+A X3 X + -/iV = o, ( 10 ) 



axvo ,xv ~dB x dA xy 2 dA xy , dA xy . . 

A xyO = ^ xy + 2 + + 2^—n 2 + -^^n y n x , 11 

oily on^ on y y on^. y 
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n x n y 



dA xy 2 
dn 2 . 7 



fl2l 



and /x x = /k x i s the background chemical potential. In evaluating these coefficients, one sets n 2 lp = n n n p after the 
partial derivatives are taken. 

In the single- fluid layers the only matter equation concerns the nuclcon radial profile n(r): 



B°n' + -iiv' 







where now 



cm z n or?/ 



(13) 



(14) 



There are several sets of "boundary" conditions that must be dealt with: at the center, at interfaces, and at the 
surface of the star. In view of Eq. requiring a non-singular behavior at the center of the star will impose that 
A(0) = 0, and consequently A'(0) and f'(0) must also vanish. This in turn implies, in view of Eq. (JTUJ), that n' x (0) 
has to vanish as well. At the surface, we will only consider configurations that satisfy n(R) = 0. A smooth joining 
of the interior spacetime to a Schwarzschild vacuum exterior at the surface of the star implies that the total mass M 
of the system is given by 



M 



f R 

An / r 2 A(r)dr 
Jo 



(15) 



and that ^(R) = 0. The metric must be continuous across the interfaces, but the matter behaviour is a bit more 
complicated [13]. We will discuss this in detail in Sec. Ill Dl 



C. The linearized field equations 



It is well-known that all non-trivial fluid pulsation modes of a nonrotating fluid star correspond to polar pertur- 
bations (often referred to as "even parity"). In the so-called Regge- Wheeler gauge [3o| . the corresponding metric 
components arc 



e v r l H (r) iujr^H^r) 

iujr^H^r) e x r l H (r) 

r l+2 K{r) 

r l+2 sm 2 0K(r) 



where Pi{9) are the Legendre polynomials. This decomposition will be applied to each layer of the star. 
The linearization of the fluid equations follows from the basic relation 

8ni = (sv + ^ x V)^< + (^ x V + ^ x V)K+ 



(16) 



(17) 



where 



A x 



dB x 

B x g^ - ^-g^g^g^pKK , 
n dB x a 

~ dn 2 9p-< 7 9i j p n x n y ' 

/ dB x dA xy 



\dn 2 



A xy — A xy ci - ci ci 



dB x 



dB y 



dA xy 



dnl Y x 



n a n p + 
dn 2 y y dn 2 



,<r„P 



(18) 
(19) 
(20) 
(21) 



These terms account for effects due to a single constituent's bulk (B x ), the presence of multiple constituents (X xy ), 
and entrainment (A xy ). This decomposition of the coefficients into these separate classes was first made by Andersson 
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and Comer Their relation to the coefficients mentioned earlier (which have been used in different applications 
0,11 HI) are 



A xy [j 



g 00 (X^oo + *4 xy oo) 



B™ = g 00 (£ x 00 + A* 



00 J 



(22) 
(23) 



Writing the nucleon four-current as = n x ii£, where g^ut^u 1 ^ = — 1, one finds that the velocity perturbation in 
the two-fluid layer is 



Sul 



where the displacement vector has components 



at ?x ' 



(24) 



6£ = e-^V-^xCrJflt 



< = -r l - 2 V*{r)-^P l e i » t 



The Lagrangian variation for each nucleon number density can be written as 

An x = 8n x + n x e- A/ V _1 W x , 
and the conservation equation Eq. |(SJ) for each particle number current yields 

An 



i -m v i Ho _ K 



(25) 



(26) 



(27) 



Thus, all matter variables can be expressed in terms of the velocity variables W x and V x . In a single fluid layer, we 
would have as the only independent matter variables the two velocity components W and V. 

The set of perturbation equations that we solve for in the superfluid layer have already been listed in (30j . but 
since our multi-fluid and multi-layer problem requires a slightly different method of solution we repeat the relevant 
equations here. First we define for each species (in analogy with Lindblom and Detweiler's [13, [H| approach to the 
one-fluid problem) the new variable 



v/2 



( x iJ + e"^V (B x n x V x + A xy n y V y 



2 

v—\)l2 : ( B *0 nxWx + A xy0 % W J 



r 



(28) 



Then we find that the Einstein and superfluid field equations yield an algebraic constraint equation 



2-1 -I 2 



-2e x - v uj 2 + e x 



I 2 + 1-2 



,.2.\ 



l-e- x 



2uo 2 1(1 + 1) x fl-e- 



8tt* 



x ) - 47rr 2 * 



Hi 
K 



+l&Tie x ~^ 2 (X n + X p ) = , 
and a system of coupled ordinary differential equations [where we use the definition Dq = B n QB P Q — (A np {j)^ 

'A' -v' l + l 



K' 



-H 



r 

e A/2 r 



Hi + —K - 1 6tt— (fi n n n V n + fi p n p V p ) 
r r 



2r 



Z + l 



A" - 8tt- 



3 A/2 



(29) 

(30) 
(31) 



-Ho + e x/2 rK 



B y 8 
n x D° 



e (A-,)/2 rXx + ^ (B x gn x ^ x + A xy °n y W y 
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A xy ° 



n x n y Dg 



,^-^ r X y + n' y (A* y °n x W x + B y °n y W y 



X x — — X x 
r 



+pt n x 
+W 2 



e v l 2 

1(1 + 1) J /2 
4 r + 2 
i/ _ J_ 

T ~ 27 



// x n x ( - - v' ) - < (B x °ti x + A xy ° % ) 
Hi 



Hn 



-(B x °n x + A x ^ % ) < 
+ i£±il e -/2^ (B x °n x K + A^° n y V y ) - e ( A ^)/ 2 — n x (£ x n x ^ x + A^n, 



y" y y> 



47re (A+,)/2M^X (^ xW ^ x + I^UyWy) + e-^' 2 



B™'n x W x 



+A*y»nyW y ) + + 



n x 

2r x r 



B x 8n x W x + A xy °n y W yj 



In the one-fluid case, the constraint equation becomes 



2-l-f 



-2e x - u u 2 + e x 



I 2 + 1-2 



8?r* 



H + 
l-er x 



2lu 2 1(1 + 1) x fl-e 
~e v 2 6 



-A 



.2 



8?r* 



8tt^ 



1- - (1- e~ A ) - 47rr 2 ^ 



A" 



+ lQne x - y/2 X = , 
and the other two equations for the metric are 



Hi' 



-Ho 



A' - v' l + l 



Hi + —K - IQir—finV 
r r 



k' = + l M^ Hl 



2r 



v' l + l 
~2 T~ 



e A/2 

K - 8tt finW 

r 



The final two fluid equations are 

e A/2 r 



w 

X' 



w, x/0 l(l + l) l + l e U-^)/2 

H + e x ' 2 rK - e x / 2 V V - — — + X , 



j 2 B 



r 2 



f 



///) 



lift 
l/_ _ 



n'B a Q n 



B° nn' 



K 



H + fin 



1(1 + 1) , ^_,_ v/a 



4 r 



Hi 



3 (A-,)/2^ Bn2M/ _ 47re (A+,)/2 W 



where 



-(A-i/)/2 



X = n 



2 P' , A '-^'^ "'"Uo 



2r 



B W 



„v/2 



-liHo + e-" /2 cu 2 (BnV) 



C"-A)/a!L (B°nW) 



(32) 



(33) 



(34) 



(35) 



(36) 



(37) 



Outside the star, the problem is reduced to solving the so-called Zcrilli equation. We refer the reader to [30( for 
details on how to match the interior and exterior solutions. 



D. Interface Layer Junction Conditions 

At the center of the star, the boundary conditions are those given in Appendix A of Comer et al [3(|, i.e. all 
functions are regular. At the surface, the conditions are the one-fluid conditions used by Detweiler and Lindblom 
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[37U38I]. The main difference here concerns the interfaces between the single-fluid and two-fluid layers. The detailed 
treatment of an interface was discussed by Andersson et al fl8j . They found that the relativistic junction conditions 
imply that the three metric perturbations Hq,Hi and K must be continuous at each interface. 

To make progress we also need to specify the behaviour of the superfluid neutron velocity, in practice W n , at each 
interface. A truly realistic model would require a detailed analysis of the phase-transition from the single fluid regime 
to the two-fluid domain. The physics of chemical equilibrium is more complicated than that of strict particle number 
density current conservation. Hence the interface behaviour is likely to be complex, with finite temperature effects 
playing a crucial role. To see that this must be the case is easy. 

The superfluid energy gap A varies with the density, cf. Fig. [TJ and for any given temperature T one would expect 
the transition to superfluidity to take place when kgT « A. In a region of "strong superfluidity", where k^T <C A, 
one can safely ignore finite temperature effects. However, as one approaches the phase transition (at points where 
ksT — » A) excitations will play an increasing role. In fact, they dictate how the two distinct fluid degrees of freedom 
couple as one transits into the single-fluid regime. This transition problem has not yet been investigated in detail. 
Hence, we must make simplifications in order to proceed. 

One can think of two natural limits. In the first, the particle reactions that drive the system towards chemical 
equilibrium are slow compared to the mode oscillation in the transition region, leading to a situation where a "su- 
perfluid" fluid element can temporarily move into the single fluid region without losing its identity. In that case, 
Wn would be freely specificable at the interface. In the opposite limit, when reactions are faster than the timescale 
of oscillation, a fluid element that moves across the interface would immediately lose its identity. This would lock 
the velocities at the interface, and as a result W n would be linked to the single-fluid velocity W. In our numerical 
calculations we assume that it is this latter scenario that applies. This is in contrast with [18|], where it was assumed 
that reactions are slow. 



III. ANALYTIC EXPANSION FOR THE LOCAL MATTER CONTENT 



The perturbation equations require as input information about the bulk, multi-constituent, and entrainmcnt effects. 
This necessarily leads to a number of coefficients to be determined. Admittedly, for those not familiar with the multi- 
fluid formalism (some of) these coefficients may seem somewhat obscure. In particular, since most of them tend to be 
ignored in discussions of the supranuclear equation of state. However, for our present purposes, the study of neutron 
star quasinormal modes, we can simplify the problem. To do this we consider an expansion of the master function 
that accounts for the fact that our background configuration is static and spherically symmetric. Alternatively, the 
expansion can be interpreted as one where the fluid velocities are small compared to the speed of light (l8l . [24| . For 
mode oscillations, this is a reasonable assumption. It is also expected to be accurate for slow-rotation models, see 
& 

Consider, for example, a term like dA xy / dn^ p . In principle, it implies that we need to expand the master function 
to 0(n^ p ). On the scale of a fluid element spacetime can be taken to be that of Minkowski. The n np term can then 
be written as 



2 1-V n -V p /C 

n nn = n n n p — ===== — ===== , (38) 



where v n and v p are, respectively, the neutron and proton three-velocities in the local Minkowski frame. When the 
individual three-velocities v ntP satisfy i> n ,p/c <C 1, then nf = n n n p up to first-order in the ratio w x /c. 



With this as our guide we write the master function in the form [la . \24 

oo 

A(nl,nl,nl p ) =J2 X *( n l> n l)( n l P - n n n pY ■ ( 39 ) 

i=0 

The bulk, multi-constituent, and entrainmcnt coefficients thus become 

oo 

A np = -Ai - £i A * Kp " Vp)'"' . (40) 

i=2 

B* = - ^ - - £ |^ K P - n n n p y , (41) 
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Bo = -^-E^Kp-^) . ( 43 ) 
<9„4 xy 9Ai ^ / \ -_ 2 

2X2-^1 (i-l)\i(nl p ~n n n p y 2 . (45) 



<9n x <9n x ^ V^ n x 

<9„4 xy 



9 < l= 3 

The *4 xy and B x coefficients on the background are 



^ xy = -Ax , B x = - — ^+ 7 ^\ 1 . (46) 

fix 



dn-y. 



This implies for the background chemical potentials that 



^ = • (47) 

In the single-fluid case, — Ao is the energy density, and so the chemical potential above is equal to that of the single-fluid 
result. We now see that the other coefficients simplify to 



d^ = d^ B *o^. (48) 
on y on x on x 



To gain a little further insight, a local, plane wave analysis [l[ indicates that the sound speeds (plural, because 
there are two) are the solutions to 

(B n c 2 s - B n °) (B p cl - B p °) + (A np c 2 s - A np °) 2 = . (49) 

We see that in addition to the expected bulk contribution to the sound speed (through B K and B x g) we have also 
entrainment entering via .4 np and the "symmetry energy" through A np g. These effects at the local level have global 
consequences. For example, Andersson et al [l8| have shown that entrainment has an impact on the oscillation mode 
spectrum by introducing so-called "avoided crossings" between the ordinary and supcrfluid modes as the entrainment is 
varied. Prix et al [35j have demonstrated that the symmetry energy leaves its imprint on slowly rotating configurations 
when the neutrons and protons are not co-rotating. 

Another dramatic effect of the entrainment is that it may facilitate a two-stream instability (39l. |40|. In principle, 
such instabilities may operate in any system with two inter-penetrating components moving at different speeds. The 
key requirement for the instability to act is that the the two fluids are coupled, and entrainment can facilitate this 
coupling. The superfluid two-stream instability is being considered as a possible trigger for glitches and appears to 
be consistent with seven years of glitch data from the fastest known young pulsar J0537 [4l| . 



IV. EQUATION OF STATE INPUT 



A. Master Function Modulo Entrainment 



In our numerical calculations we will compare two relatively simple paramctrisations for the equation of state. 
Nothing prevents us from considering other perhaps more realistic models, including tabulated equation of state data, 
but at this point we are more concerned with the technical developments than with the actual numerical results. As 
we will discuss below, see Sec. IVI1 a fully consistent model requires information that is not yet available. Until this 
changes, this kind of analysis must be viewed as somewhat qualitative. 

We consider two models for \ n : i) a simple "TOY" model that is a sum of polytropes for the neutrons and protons, 
a symmetry energy term, and a final term that accounts for a relativistic gas of electrons, and ii) the more realistic, 
so-called "PAL" equation of state [42j ■ 
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Model 


T (MeV) 


n(0) (fnT 3 ) 


M (M s ) 


R (km) 


TOY 




1.0 


1.635 


10.655 




0.45 


1.0 


1.745 


10.641 


PAL 


0.5 


1.0 


1.751 


10.661 




0.55 


1.0 


1.759 


10.684 



TABLE I: The "canonical" background stellar models for the TOY and PAL equations of state for the parameter values 
discussed in the main text. 



Specifically, the TOY model in the one-fluid and two-fluid layers takes the form 

A /m = - (n„ + n p ) - o h (n n + n p f b - S (n n + n p ) (1 - 2x p ) 2 - a e n^, B , 



(50) 



where we have chosen the parameters at, = 0.2, a c = 0.5, /3b = 2.5, j3 c = 2.0, So = 0.05, m is the baryon mass, 
and x p = n p /n is the proton fraction. As shown in Table [H this combination leads to a neutron star model with 
reasonable mass and radius. 

The particular PAL master function for both layers is 



Xo/m = — tiqu 1 + Eq(u) + S(u) (1 — 2x P Y 



A /m , 



(51) 



where 



mE (u) = A u 2/3 + Bou + Cqu* + 3^2^^ 
mS{u) = A s (u 2/3 - uj + S u 



,=i 



,1/3 



tan 



1 (a^ 1 



/3 



A, 



— X(Xe ) 
T 



Stt 



^ {x (1 + x 2 ) 1/2 (1 + 2x 2 ) - In [x + (1 + x 2 ) 1/2 ] } , 

r e = h/m e , 



1836 



3tt^ 



1/3 



,1/3 



(52) 



with u = n/n (n = 0.16 fm" J ), u = 0.927, A = 22.11 MeV, B = 220.47 MeV, C = -213.41 MeV, C x = 
-83.84 MeV, C 2 = 23.0 MeV, cti = 2/3, a 2 = 1/3, A s = 12.99 MeV, and S = 30 MeV. Again, the values for mass 
and radius obtained for these parameter values are reasonable, sec Tabic |U 

The electron contribution A is vital for ensuring chemical equilibrium of the system. Even though the electron 
mass m e is much smaller than the nucleon mass (m e = m/1836), the fact that the electrons are ultra-relativistic gives 
them enough energy to affect the chemical potential at the same level as the nucleons. Because of overall charge 
neutrality we set n c = n p , and recall that the electrons and the protons flow together. Note also that, even though the 
imposition of chemical equilibrium will effectively make the single-fluid and two-fluid equations of state the same on 
the background, we have to distinguish them because of the £> x , B x g, etc. coefficients, for which the partial derivatives 
must be computed before chemical equilibrium is imposed. 



B. Gap Model 



BCS theory is the basic paradigm for Fcrmionic superfluidity Given a many-particle system with an attractive 
component in the interactions, the theory tells us that Cooper pairs will form, leading to a fundamental modification 
of the energy states near the Fermi surface, namely the formation of a gap in the energy spectrum. It is the existence 
of this gap that leads to superfluidity and in the present context multi- fluid dynamics. Even if particles try to scatter 
dissipatively, there are essentially no accessible states for them to scatter into unless there is enough energy to break 
a Cooper pair. When pair-breaking occurs the particles can reach the states above the gap, and energy can be 
irreversibly dissipated. Those particles in energy states "below" the Fermi surface are in this sense also "supcrfiuid" 
since they cannot enter already filled states. So, for example, when neutrons are superfluid, ordinary scattering that 
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would cause (on average) neutrons and protons to flow as a single fluid is severely diminished. It is thus clear that 
neutron star dynamics on a macroscopic scale is dictated by the gap structure on the microscopic scale. 

There has been much work analyzing the details of the gap structure in dense, nucleonic matter trying to account 
for different interactions, medium effects etcetera (see @ for a review of gaps in general and Q for an application 
to neutron stars). If even our rudimentary understanding of the strong force is correct for the densities expected in 
neutron stars there is little doubt that neutron and proton gaps exist. But given the complexities of the problem there 
is no exact agreement on the details; particularly on the density dependence of the gap. This affects the maximum 
gap energy, the gap profile as a function of density, and the size of regions in which gaps exist. However, it is generally 
accepted that neutron superfluidity and proton superconductivity will not extend throughout the star. Consequently, 
there should be layers having ordinary, single-fluid dynamics as well as ones where the multi- fluid description applies. 
This implies that most, if not all, previous mode calculations and modeling of rotational equilibria that assume 
superfluidity throughout the star are incomplete (see [l|, HH for reviews). 

As a first step towards improving the situation we will consider models that account for the detailed energy gap 
and its dependence on, in particular, the density. For practical reasons, we will use the parametrized model of 
piil . 0, [ill, SH| as adapted by Andersson et al [g]. This model has a number of "free" parameters, that can be 
used, for instance, to adjust the maximum gap energy Ao and the gap profile (as a function of density). In this 
phcnomcnological description, the gap energy A(fcp) (where /cf is the wave- number at the Fermi surface) takes the 
form 

A(M = A (fcF ~ fc 2 l)2 y 2 ■ (53) 

(fci? - k\) + k 2 (&f - fci) + ki 

Note that density dependence enters implicitly through the wave-number kp. Andersson et al [f| provide a table of 
parameter values representative of the various gap calculations in the extant literature. For our numerical calculations 
we focus on gap model "h" given in their Table 1. 

Once an equilibrium configuration is built, and density is known as a function of radius, we can use the gap energy 
to determine which layers of the equilibrium configuration are ordinary or superfluid. This requires an assumption 
of the temperature profile in the star. In principle, our calculations assume that the fluid is at zero temperature, 
somewhat in the same spirit as mode-damping calculations due to shear and/or bulk viscosity [47j]. We certainly 
do not have a consistent temperature description. However, for our present purposes, this is not a major problem. 
The gap information that we use is phcnomcnological so it should be acceptable to account for the temperature in an 
approximate way as well. In view of this, we simply assume that the fluid is isothermal (and do not account for the 
gravitational redshift). 

It would certainly be possible to improve on this description and it will eventually be important to do so. However, 
one would then like to account for all temperature effects, including quasiparticle excitations in the superfluid. This 
is an interesting problem since it opens the door for studies of dissipative superfluid neutron stars, but we will not 
discuss it further here. 

Once we have chosen a core temperature we can work out if there are superfluid regions in the star. Varying the 
temperature thus affects the size of the superfluid layer, but has only a small effect on global parameters like the mass 
and radius, see Table [T] We also adapt the number of independent fluids for the mode calculations accordingly. For 
our canonical TOY and PAL neutron star models, the particle number density and gap energy as a function of radius 
typically take the form shown in Fig. [TJ As anticipated, there is a layering of regions having different fluid dynamics. 
Wherever the gap energy A is greater than the thermal energy ksT, the neutrons will be superfluid, and therefore 
two-fluid dynamics will apply. 



C. The a - co Relativistic Mean Field Model 



The model to be used for entrainment has been obtained using a relativistic a — tu mean field model of the type 
described by Glendenning [ll[ ■ Whatever reservations one may have about mean-field equations of state, the great 
advantage from our present perspective is that the entrainment can be quantified. The Lagrangian for this system is 
given by 

L = L b + L a + L u + L int , (54) 

where 



Lb = ^(il^ - m)ip , 



(55) 
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FIG. 1: Plot of the TOY model density profiles (vs radius) and gap function (vs the Fermi momentum) for the case ksT = 0.4 
MeV. The dashed line is the superfluid neutron density, and in the same range the solid line is that of the protons. Otherwise, 
the solid line represents the total baryon number density in the non-superfluid regions. If the gap energy is higher than the 
thermal energy, the neutrons are taken to be superfluid. 



La 

Lint 



1 ui/ 1 2 u 

— u)„ v u)^ m,,(jj,,w 



(56) 

(57) 
(58) 



Here m is the baryon mass, ip is an 8-component spinor with the proton components as the top 4 and the neutron 
components as the bottom 4, the 7 M are the corresponding 8 x 8 block diagonal Dirac matrices, andw M „ = d^uj^ — d^uj^. 
The coupled set of field equations obtained from this Lagrangian are to be solved in each fluid element. 

The main approximations of the mean field approach are to assume that the nucleons can be represented as plane- 
wave states and that all gradients of the a and fields can be ignored. The coupling constants g a and g^ and field 
masses m a and are determined, for instance, from properties of nuclear matter at the nuclear saturation density. 
Fortunately, in what follows, we only need to provide the ratios (? a = {ga/ma) 2 and c£ = (g^/m^) 2 - 

The main aim here is to produce a master function that incorporates the entrainmcnt effect. Consider again fluid 
elements somewhere in the neutron star. The fermionic nature of the nucleons means that they are to be placed into 
the various energy levels (obtained from the mean field calculation) until their respective (local) Fermi spheres are 
filled. The Fermi spheres are surfaces in momentum space. The entrainmcnt is incorporated by displacing the center 
of the proton Fermi sphere from that of the neutron Fermi sphere. The neutron sphere is centered on the origin, and 
has a radius fc n . Displaced an amount K from the origin is the center of the proton sphere, which has a radius fc p . 
The Fermi sphere radii and displacement (k n , k p ,K) are functions of the local neutron and proton number densities 
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OJ 0.4 



FIG. 2: The radial profile of e mom defined in Eq. (for <* = 7.148 and cl = 12.684) using the same TOY model as in Fig. ED 



and the local, relative velocity of the protons, say, with respect to the neutrons. 

Earlier we discussed the analytic expansion that applies to the master function, and determined that we need to 
know the Aq and Ai coefficients. The mean field model yields [24 1 



Ar 



■5// 



ii 2 



2A- 



2 \J fc 2 + m* 



3tt 2 



h 2 h 3 



3^ 2 fc 2 k 2 n + m* 2 



where the Dirac effective mass m* is the solution to the transcendental equation 

,,2 



m — m* 



2tt 2 



k n \/k 2 + m* 2 + fc p J k 2 + m» 2 — 



— TO* 2 ln 



fc n + k\ + m t 



-k n + ^/k 2 + m* 2 I 2 



— m* 2 ln 



fcp + J fcp 



fcp + wfcp 



and each nucleon number density n x is related to its Fermi surface radius fc x via 

n x = ■ 
Sir 2 

Typical behaviour of the cntrainmcnt parameter e mom used in several studies, defined as 

m A np 



,2 ' 



(59) 



(60) 



(61) 



(62) 



is shown in Fig. O We set c 2 = 7.148 and c 2 = 12.684 in what follows. These values are consistent with known 
nuclear matter properties and generate reasonable neutron star models [ll|, H3 ■ 



V. NUMERICAL RESULTS 



We consider three variations on the mode calculations: i) changing the master function, ii) including (or not) 
cntrainmcnt, and iii) varying the temperature. As we have already discussed, the main effect of the latter is to alter 
the size of the two-fluid layer. Since the parameter space is vast, and we are mainly interested in understanding 
the qualitative effects, we restrict our discussion to a few "canonical" background models. These determine the 
total mass M and radius R given a particular central baryon number density n(0) = ?^ n (0) + ?^ P (0), obtained by first 
specifying rt n (0) and then finding n p (0) through the condition of chemical equilibrium. Of course, as we vary the 
temperature T of the star, the matter distribution is affected which, in principle, changes the values of M and R. 
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FIG. 3: Plot of log 10 |A in | ve Re{u;M) for the TOY model, using k B T = 0.38 MeV and e mom = and s mom + 0. 
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FIG. 4: Plot of log 10 |An| vs Re(wM) for the TOY model, using k B T = {0.38,0.4,0.42} (MeV). Note the "gap" modes 
(indicated by arrows) that shift as the temperature (i.e. the thickness of the two-fluid layer) is varied. 



However, we find that these changes are negligible for the TOY model, and only slight for the PAL model (cf. Table U 
and the disscusion in Sec. lIVBj) . 

A quasinormal mode corresponds to a purely outgoing (gravitational) wave solution to the perturbation equations. 
The asymptotic amplitude for ingoing waves, Ai n , can therefore be used to locate the modes. Figs. [3J QJ El and [5] 
show this amplitude versus frequency uiM for different models (see f30l ] for the definition of Aj n ). The real part of 
the different quasinormal-mode frequencies can be approximated by the values that make Ai n vanish (or log 10 l-<4-m| 
tend to -co). This technique does not provide the damping times of the modes due to gravitational- wave emission. 
In practice, we can use the technique described by Andersson et al [l|| to reliably determine the damping times. Our 
present discussion will, however, focus on the oscillation frequencies. Figs. [Q [H and [9] provide the radial profiles of 
mode cigenfunctions for the first few frequencies. 

Several features of the mode spectrum arc important to note. As in previous studies [Til . |30T ] of the general 
relativistic two-fluid system, there are "more modes" than in the single-fluid case. This is expected since we have an 
additional fluid degree of freedom (albeit localised to a distinct layer in the star). From Fig. [3] we see that entrainment 
has a small effect on the first mode in the spectrum, the f-mode, but clearly affects the next, a "superfluid" mode, 
frequency by shifting it from the value labelled u\ to that labelled v\. Such behaviour is consonant with the local 
analysis of Andersson and Comer (48j . who show that the additional, superfluid modes depend on entrainment, while 
the fundamental f-mode does not. 
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FIG. 5: Plot of log]A in j vs Re(wM) for the PAL model, using k B T — 0.5 MeV and for e mom = and e mom 7^ 0. 
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FIG. 6: Plot of log|A m | vs Re(wM) for k B T = {0.45, 0.5, 0.55} (MeV) for the PAL model. 

From Fig. [4] we see that as the size of the two-fluid layer is decreased most of the mode frequencies remain 
largely unchanged. However, in this model there are clear "gap" modes (emphasised by the arrows) whose defining 
characteristic is that the layer thickness (i.e. the size of the gap energy A relative to the temperature ksT) has 
significant impact on the frequency. As might be expected the frequencies increase as the thickness decreases. This 
is natural since the mode wavelengths must adjust to the size of the "container" . 

For a semi- realistic determination of mode frequencies we consider Figs. [5] and [5] which show fog 10 |A;„| vs loM for 
the PAL master function. Qualitatively, Fig. [5] shows that the effect of entrainment in the PAL model is the same as 
for the TOY model. Perhaps more interesting is the observation that the layer thickness affects all modes in the PAL 
model, cf. Fig. [6l However, the most pronounced change is that the second "gap" mode (indicated by the arrow on 
the right of the figure) occurs earlier in the mode spectrum than in the case of the TOY model. Hence it seems clear 
that the combined background master function and the gap data are to some extent distinguishable via the mode 
spectrum. 

As expected, the radial mode profiles have a richer structure than in the single-fluid case. There are also marked 
differences with models where the two fluids extend throughout the core [l8|, H(| ■ In each of Figs. [7J [5J and[51 the solid 
line represents the total baryon flow in the one-fluid regions. In each two-fluid region the dot-dashed line represents 
the neutron flow, the dotted line represents the proton flow, and the dashed line represents the net baryon flow. 
Note that the boundary conditions between the one- and two-fluid regions imply that the net baryon flow must be 
continuous. 

In Figs. [7J and [5] we see a characteristic feature of superfluid modes, the counter- motion between the neutrons and 
protons. The final Fig. [9] illustrates the effect of entrainment. In the plot for the vs mode we see the beginnings of 
an "avoided" crossing. As first demonstrated by Andersson et al [l8| , when viewed as a function of the entrainment 
two neighbouring mode frequencies sometimes come close to crossing but veer off before intersection happens. Past 
this avoided crossing the two modes exchange character in the sense that a co-moving mode becomes counter-moving 
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FIG. 7: Plot of the TOY model eigenfunction radial amplitude (vs r/R) for the first "gap" mode (i.e. the arrow on the left in 
Fig. [4} for the case UbT — 0.38 MeV. The frequency of the mode is Re(cjAf) = 0.3583. The mode is such that the two fluids 
exhibit counter-motion in the radial direction; i.e. as the neutrons move in, say, the protons are moving out. 
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FIG. 8: Plot of the TOY model radial amplitude for the second "gap" mode (i.e. the arrow on the right in Fig.|4]) for fcsT = 0.38 
MeV. The frequency of the mode is Re(wM) = 0.6066. The mode is such that the two fluids exhibit counter-motion in the 
radial direction; i.e. as the neutrons move in, say, the protons are moving out. 

and vice versa. 



VI. DISCUSSION: NEXT GENERATION EQUATIONS OF STATE 

In this work we have extended the general relativistic two-fluid formalism used to model quasinormal-mode oscil- 
lations in such a way that superfluid neutron gap data can be incorporated. The main point is that the presence 
of a gap not only increases the number of fluid degrees of freedom locally; it also introduces a dynamical layering in 
the sense that the number of fluid degrees of freedom can change depending on whether the thermal energy of the 
star in a region is greater or less than the gap energy. It is clear that this structure will also be important for stellar 
rotational equilibria (eg. as determined by a slow-rotation approximation). 

In our models we analyzed the effects of different master functions, entrainment (obtained from a relativistic a — lo 
mean field model), symmetry energy, and the presence of a gap. We found that the effect of entrainment on the 
mode spectrum is qualitatively consistent with earlier studies. On the other hand, the combined effects of the gap 
and the master function have distinguishable effects on the mode spectrum. Basically, our results illustrate that any 
"realistic" model for oscillating superfluid neutron stars must account not only for the "bulk" equation of state, one 
must also incorporate the entrainment and the superfluid gap data in a consistent way. We have, of course, not 
done this. At this point our models arc more or less phenomcnological. However, we have demonstrated that the 
computational technology required for this kind of study is now in place. What is missing is the input microphysics. 
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FIG. 9: Plot of the TOY model eigenfunction radial amplitudes (vs r/R) for the modes «i..,ui...etc labelled in Fig. EI The 
dashed, dot, and dot-dashed lines have the same denotations as given in the insets of the previous two figures. The remaining, 
horizontal line is to pinpoint where the amplitudes vanish. 



At the present time there is, as far as we are aware, no single equation of state calculation that provides all the 
required parameters. This needs to change in the future. Furthermore, one must not forget the neutron star crust, 
where the lattice of nuclei will have elastic properties. These also need to be determined consistently. 

Let us discuss what (in our view) is required from the next generation equations of state if we want to model 
mode oscillations of multi-fluid neutron stars. To some extent these demands should be relatively straightforward to 
meet. Like most complex models, realistic equations of state are presented in tabular form. The question is what 
to include in the table. For traditional reasons, the information usually given is the pressure as a function of density 
and/or temperature. This would be adequate if one were only interested in solving for a background configuration 
single-fluid star. The information is, however, not even complete for the oscillations of such models. In order to 
model the so-called g-modes, which are sensitive to both temperature and composition gradients, one would need also 
the various particle fractions. Moreover, the multi-fluid formalism presented here requires a number of additional 
quantities. 

In our experience, the two-fluid problem is most simply understood as having two independent thermodynamic 
parameters, the two number densitites {n n , n p }. We see from the background and perturbation equations in Section HI1 
that we require the set of values {A, 'J, B K , A xy , d^ x /dn y , A x }. Each variable is a function of {n n ,n p }, and hence 
they will need to be given as two-dimensional arrays, the first index of each array standing for n n , say, and the second 
for n p . Of course, one usually assumes chemical equilibrium among the particle species on the background, which 
implies, say, that n p = n p (n n ). However, a number of the quantities that are needed for the perturbation analysis, 
like the entrainment A xy , require information about the system away from equilibrium. If we let T represent the set 
of thermodynamic variables, then a sample tabular equation of state could look something like that of Table ITT1 

Can this information be provided within a single framework for determining the "equation of state" ? We do not 
see why it should not be possible. Some of the quantities wc require arc already calculated, they arc simply not 
presented in the final equation of state table. It should certainly be very easy to include the particle fraction in the 
tabulated data, and we see no reason why the entrainment coefficients should not be straightforward to determine 
as well. It is obviously the case that there is no universally agreed upon equation of state, or indeed method for its 
determination. From our point of view this is less relevant. We expect to live with uncertainties. After all, we 
do not have the luxury of having neutron stars readily available in the laboratory. The key point is that consistent 
models require all parameters to be consistent with a given microphysics calculation. Hopefully, an improved dialogue 
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TABLE II: A schematic representation of a realistic, tabular equation of state that gives N x N entries for the set of thermo- 
dynamic variables T = {A, *, S x , A* y , dfx x /dn y , A x }. 

between neighbouring areas of research will allow us to make progress. 

There are, of course, a number of related challenges. We have not yet satisfactorily answered the question of 
the minimal number of independent fluid degrees of freedom required to model a compact star. Glitch data tells 
us already that there are at least two. But what if we think about the deep core, quark dcconfincmcnt and CFL 
matter [lH, HH ? Do the various Cooper pairings between quarks indicate the presence of independent "fluids" ? If 
they do, then what physical mechanisms exist to excite the additional degrees of freedom, and what is their physical 
interpretation? How can one hope to calculate the myriad of phases that could occur and then translate that into 
numerical models of rotating and oscillating multi-fluid compact stars? These are challenging questions that require 
further consideration. 
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Appendix: Boundary conditions and numerical scheme 

We have three different regions inside the star: (1) one-fluid core; (2) two-fluid superfluid region; and (3) one-fluid 
envelope. In the one-fluid regions, fi is the baryon number density (with /2 being the chemical potential). In the 
superfluid region, n n and n p are the neutron and proton densities respectively (n = n n + n p ). Note that, in this 
Appendix, we use tildes to make a distinction between the variables in the single-fluid and the multi-flud regime. Ri 
and i? 2 denote the corc-supcrfluid and superfluid-crust interfaces, respectively. We let R denote the star's radius. 

The one-fluid regions are governed by 4 ODEs for 

Y = {H 1 ,K,W,X} . (63) 
The superfluid region is governed by 6 ODEs for 

Y = {H U K, W n , W p , X n , X p } . (64) 

A. Boundary conditions 

At r = 0, we have two conditions relating {i?i(0), K(0), W(0), X{0)} (cf Eqs. (A5)-(A6) of Comer ct al [H, with 
n a = 0, Wn(0) = 0, n p — > n, etc. in our current notation). Explicitly, we have 

„f /2 / ,,2 \ 

X(0) = —^-jxnK{0) - ( e Uo/2 n 2 B°h + —e- Uo/2 Bn 2 J W(0) (65) 

and 

ffx(0) = j^jK(0) + -^mW{Q) . (66) 
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At each of the interfaces (i?i and R2), we have the following 4 junction conditions according to (A16)-(A18) of 
Andersson et al [T^j: 

Hi{R c ) = H X {R C ) , (67) 
K(R C ) = K(R C ) , (68) 
fi(R c )n{R c )W(R c ) = fi n {R c )n n {R c )W n {R c ) + n P {R c )n p (R c )W p {R c ) , (69) 
X(Rc) = X n (R c ) + X P (R C ) , (70) 

where we have assumed that A and VP are continuous across the interfaces. 

We further impose that the two fluids move in lock-step at the interfaces. This effectively translates to the condition 

W n {R c ) = W P (R C ) . (71) 

In summary, Eqs. (f6"7|) - (f7Tj) are the required boundary conditions at the interfaces. 
Finally, at the surface of the star r = R, we have the single condition 

X{R) = . (72) 



B. Numerical scheme 

At r = 0, we need only to specify {K(0),W(0)}. The remaining variables {-ffi(O), X(0)} are determined by 
Eqs. (l65l) - (|66|) . All of the second derivatives Hi (0), K"(0), etc. are also determined. We choose two arbitrary 
values of {K(0),W(0)} and integrate the 4 ODEs from small r up to the core-superfluid interface r — R\. The 
general solution in the core region is 

2 

Y(r) = f,Y,(r) , for < r < R 1 . (73) 

i=X 

Next we turn to the general solution in the superfluid region R\ < r < i? 2 ■ At r = Ri, the solution must satisfy 
the condition (|7Tj) . This means that we must generate 5 linearly independent solutions. This is obtained by choosing 
five different sets of {-ffi(-Ri), K (R{), W n (R x ), W p (Ri), X n (R{), X p (Ri)} (with W n (Ri) = W P {R X )) and integrating 
to r = i?2- The general solution in this domain is 

7 

Y(r) = CiYi(r) , for R t < r < R 2 . (74) 

i=3 

At the surface, our solution must satisfy X(R) = 0. Hence, we must generate 3 linearly independent solutions in 
the crust (i?2 < r < R). The general solution in this domain is 

10 

Y(r) = Y CiYi(r) , for R 2 < r < R . (75) 

i=8 

After fixing the overall normalization by choosing the value of one of the c,; (says, C10), we have 9 remaining constants 
Cj (i = 1, .., 9) to be determined by the boundary conditions at R\ and R 2 . 

First, at r = we have 4 conditions Eqs. (f6"7} - (|70[) to be satisfied. Note that the condition (jTTj) has been used 
to generate the general solution in the superfluid region. Explicitly, these conditions become 



5>iTi W k = 5>2?i (<) k , (76) 

i=l i=3 
2 7 

Y,CiK^\ Rl = J2*K (i) \ Rl , (77) 

i— 1 i— 3 

2 7 

Y^ (jinW^) \ Rl = (w*W® +M P n p W^) | Bl , (78) 

?—l 1—3 
i— 1 ? — 3 
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where we have defined Y. t = {H^ , K® , w£\ , X%\x^} and Yj = {Hi®, K®, W^, 
Next, at r = R2, we have 5 conditions Eqs. ([67| - (I7T1) : 



^CifwP-W®)^ = 0, (80) 



i=3 



7 10 

= E c -^ w k . ( 81 ) 



i=3 i=8 
7 10 



E** (<) k = E c * /?w k . ( 82 ) 

i=3 i=8 
7 10 

E c 4 (A*nWi i5 + lipnpWW) k = E c * k , (83) 

i=3 i=8 
7 10 

E Ci + 4°) k = E C ^ W k ■ (84) 



i=3 i=8 



The 9 equations Eqs. |76|) -([84 |) can be used to determine the 9 constants c, (i = 1, ..,9). This completes the interior 
problem. 
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